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We construct a rejection-free Monte Carlo method for the hard-disk system. Rejection-free Monte 
Carlo methods preserve the time-evolution behavior of the standard Monte Carlo method, and 
this relationship is confirmed for our method by observing nonequilibrium relaxation of a bond- 
orientational order parameter. The rejection-free method gives a greater computational efficiency 
than the standard method at high densities. The rejection free method is implemented in a shrewd 
manner using optimization methods to calculate a rejection probability and to update the system. 
This method should allow an efficient study of the dynamics of two-dimensional solids at high 
density. 
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I. INTRODUCTION 

Monte Carlo (MC) methods have become more pow- 
erful tools with the development of faster and more ac- 
cessible computers. Many different phenomena have been 
studied with MC methods _,li.2j. The melting behavior of 
the hard-disk system is one of such subjects which has 
been studied with Monte Carlo methods H S 0, S H, . 
Consider the hard-disk system with TV particles. The 
Monte Carlo method for the hard-disk system has the 
following steps. First, choose one particle randomly, i.e. 
choose it with a probability of Then choose a new 

position for the center of the chosen particle, the new 
position is chosen uniformly in the circle with radius as 
centered on the original position of the particle. The trial 
move is accepted when the new position has no overlap 
with other particles, otherwise, the trial move is rejected. 

The MC methods introduced above have been used to 
obtain the equilibrium state of the system. To study the 
equilibrium state of the system, the value of as, the trial 
step length, is often chosen in order that the rejection rate 
is near 50%. Recently, Watanabe et al. [wL fill studied 
the hard-disk system by observing nonequilibrium behav- 
ior using a molecular dynamics (MD) method. In order to 
compare the results of MD and MC, we have to consider 
the dynamics of the MC method. The dynamics of the 
MC method for the hard-disk system can be understood 
as a Brownian motion, like pollen particles in a fluid. The 
value of as corresponds to the amplitude of the external 
random force. Therefore, it determines a time scale of 
this system. If values of as are not identical for all den- 
sities, we cannot compare the dynamics, like relaxation 
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phenomena and fluctuations, between different densities. 
However, keeping as fixed leads to high rejection rates in 
the Monte Carlo at high densities. A similar problem oc- 
curs for spin systems both at low temperatures and in a 
strong external field: a rejection rate becomes very high, 
so a huge number of trials is required to make a change. 
This problem can be avoided using a different technique 
called the rejection- free Monte Carlo (RFMC) method. 
The RFMC method was first constructed for discrete spin 
systems |l2l j , for a review see jl^] . Miyashita and Takano 
[Tjl applied the RFMC method to the kinetic Ising model 
in order to s tudy dynamical critical behavior. Recently, 
Mufioz et al. [T^ proposed a new RFMC algorithm which 
can treat models with continuous degrees of freedom. In 
this paper, we construct and utilize a RFMC algorithm 
for the hard-disk system based on the method of Ref. . 



II. METHODOLOGY 
A. A Rejection-free Monte Carlo method 




P(l|2) 

FIG. 1: A Markov chain of Monte Carlo steps. 



2 



A Monte Carlo method is an implementation of a 
Markov process on a computer, and hence is sometimes 
called a Markov Chain Monte Carlo. The Monte Carlo 
method calculates various physical quantities by updat- 
ing states of a system using random variables. These 
updating processes can be illustrated by a Markov chain 
(see the schematic in Fig. ^ . Let the current state be 
at Sq and the states possible to move from Sq are de- 
noted hy Si{i = 1,2, Ns). Define Ei as the energy of 
the state Si. The new state Si{i = 0, 1, • • • , N^) will be 
chosen with probability P{i\0). One way to ensure that 
the system will relax to the equilibrium state is to insist 
that the probability P{i\0) satisfy the detailed balance 
condition 0. 

One of the well-known ways ^ to satisfy the detailed 
balance condition is to use a heat-bath transition prob- 
ability. In the heat-bath method, the probability P{i\Q) 
is defined to be 



exp i-pE,) 
Y.k=o exp{-(3Ek) 



(1) 



. When a system has a continuous degree of freedom, 
the summation of Eq. becomes an integration which 
is generally difficult to calculate analytically. 

Another popular way to satisfy the detailed balance 
condition is a Metropolis method. In this method, each 
step contains two parts; selecting a new state and accept- 
ing or rejecting the trial to move to the selected state. 
First, pick up a state Si from all possible states to move 
to with uniform probability l/Ns- The probability P{i\0) 
to move from 5*0 to Si is defined to be 1 when AEi < 0, 
otherwise it is exp {—f3AEi) with the energy difference 
AEt = Et - Eo- Therefore, the probabihty P{i\0) in the 
Metropolis method is 



P(z|0) - 



l/Ns 
exp{~pAE,)/N, 



if AE,<0, 
otherwise. ^ ' 



When a trial is rejected, the configuration of the system 
is not updated. The probability A = f'(OlO) to stay in 
the current state after the trial is given by the expression 



A = 1 



^P(^IO). 



(3) 



For some parameters, e.g. under a strong external field 
and at an extremely low temperature, the value of A can 
be very nearly 1. In such cases, most of computational 
time is spent on calculating trials which will be rejected. 
This rejection rate drastically decreases the efficiency of 
the computation. 

In order to overcome this problem, a rejection-free 
Monte Carlo (RFMC) method is proposed. It is an ex- 
ample of an event driven algorithm [ll, l2l and has also 
been called a waiting time method [Tj, • Each steps 
of the RFMC method involves first computing the time 
to leave the current state (the waiting time t.^g^n), and 
then choosing a new state to move to with the appro- 
priate probability. It does not contain the judgment to 



accept or reject a trial, and, therefore, it achieves rejec- 
tionless updates of the system in each algorithmic step. 
The waiting time i^^^ait calculated using a (pseudo) 
random variable. The probability p{t) to keep staying at 
the current state for t steps decays exponentially. 



p(t) = A* = exp (tin A), 



(4) 



with A defined in Eq. © . Note that in A < since < 
A < 1. Inversely, the time t to stay in the current state 
is determined with a uniform random number r on (0, 1) 
to be 



t. 



wait 



Inr 



I, 



(5) 



where [x\ denotes the integer part of x and the rounding 
down is introduced to express the discrete time step in 

Mc [mii. 

After the time of the system is advanced by t-^^ait' ^ 
new state Si is chosen from the all states possible to 
move to except for the current state with the probabil- 
ity proportional to P{i\0) 0,0,0. Since all values 
of P{i\0) are required to proceed one algorithm step in 
the RFMC, the computational cost of one step is higher 
than the normal MC. However, the waiting time t.^g^[^, 
the time which can be advanced in one algorithmic step, 
becomes large at low temperatures, and consequently the 
efficiency of the RFMC can become better than that of 
standard MC. 



B. Application to hard-disk systems 

Consider a hard disk system with N particles with the 
radius a. A standard MC method for the system involves 
choosing a particle, and trying to move the chosen par- 
ticle within a circle with radius CTs centered on the orig- 
inal position of the chosen particle. To apply a RFMC 
method to the hard-disk system, define Xi as the prob- 
ability that a trial to move particle i is rejected (given 
that particle i was chosen as the particle to attempt a 
move). Using the definition of A^, we can construct the 
algorithm of the RFMC method for the hard-disk system 
as follows: 

1. Calculate the waiting time t^^it using Eq. Q with 

2. Advance the time of the system by i-^^f^it' 

3. Choose a particle i with the probability propor- 
tional to 1 — Xi, which is the probability that (given 
that particle i was the particle chosen for an at- 
tempted move) the trial to move the particle i 
would be accepted. 

4. Choose the new position of the chosen particle i 
uniformly from all the points to which the particle 
i is allowed to move. 
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The steps described above are the same as the RFMC for 
continuous spin systems 0, but the algorithms to cal- 
culate Xi and to determine a new position of the chosen 
particle are unique to the hard-disk system. The proba- 
bility Xi can be calculated to be, 



with an area Ai into which the particle i is allow to move 
(see Fig. Therefore, if we can calculate Ai for all 
particles, then we can construct the RFMC algorithm 
for the hard-disk system. 




FIG. 2: (Color online) A schematic drawing of the definition 
of Ai (shaded). The solid circles are particles and the small 
dashed circle has a radius as- The shaded area is the area 
which is a continuous set of the points that the center of the 
chosen particle can move to. The ratio of Ai to the area of the 
trial circle ivas^ gives the probability of accepting the move, 
1 — Ai, given that the center particle has been chosen as the 
one to move. 



C. Calculation of Ai 

The area Ai is the continuous set of positions in which 
the particle can be placed without any overlaps. Without 
neighboring particles, the shape of Ai would be a filled 
circle with a radius CTs- Let's call it a trial circle. In the 
general case, the shape of Ai is the remaining part of 
the trial circle after removing the overlap of 'shadows' of 
neighboring particles. The shape of the shadow is a circle 
with a radius 2a which is concentric to a neighboring 
particle. Let's call this a shadow circle. The area Ai, 
thus, consists of areas of arcs of a trial circle and that of 
shadow circles. 

To compute the value of Ai, we develop a method we 
call the survival point method. See Fig.|31(a). The cho- 
sen particle is shown as a solid circle, the trial circle is 
shown as a concentric dashed circle, and the area Ai is the 



shaded region. Each neighboring particle (filled circles) 
has a shadow circle which is concentric and has radius 2a. 
An enlargement of the area Ai is shown in Fig. O (b). It 
is seen that in this example this area has five vertices 
which are intersection points of shadow particles, we call 
them survived vertex points. In Fig.|31(c), these survived 
vertex points are shown as small filled circles. Straight 
lines connect the center of the chosen particle and the 
intersection points. In this case the area Ai is divided 
into five figures. 

Each divided figure is the remaining part of an isosce- 
les triangle with the overlap of a shadow circle removed. 
It is easy to calculate this area. Thus, all we have to do 
is to find all survived vertex points which form the area 
Ai. First, make a list of all intersection points of shadow 
circles and the trial circle. Next, remove points which 
are included in other shadow circles from the list, since 
these points cannot be vertices forming the area Ai. Af- 
ter this removal process, we have the vertices which form 
the area Ai (see Fig. |31(c)). The calculation process of 
a partial figure which forms Ai is shown in Fig. O (d) . 
The vertices are denoted by Pi and P2, and the center 
of the shadow circle is denoted by S. The survived ver- 
tices Pi and P2 are on the shadow circle centered at S, 
so ST\ = Jl^ = 2(7. The area of OP1SP2 can be calcu- 
lated by summing the two triangles OP1P2 and SP1P2 
with Heron's formula. The area of the chord is Aa^O. 
Finally the portion of the area OP1P2 is calculated by 
subtracting the area of the chord SP1P2 from the area of 
the quadrilateral OP1SP2. The total area Ai is the sum 
of one such calculation for each survived vertex. 



D. Choosing a particle to move 

After calculation of t.^g^[i and advancing the time of 
the system by it, we have to choose a particle i to move 
with a probability proportional to 1 — Ai. With a direct 
implementation, i.e., with the integration scheme |l8j| . 
the order of the computation is 0{N), which is very time 
consuming. Other approaches are proposed like a three 
level search for spin systems ,19|. The three- level search 
improves the efficiency of the search by determining co- 
ordinates of a spin to update one by one. However, it is 
difficult to apply this method for particle systems, since 
neighbors of particles are not fixed. Here we use a com- 
plete binary tree search for the choosing part of the al- 
gorithm. 

First, calculate the area Ai for each of the particles. 
Since an acceptance probability 1 — Ai is proportional to 
Ai as shown in Eq. ®, the particle should be chosen 
with the probability proportional to Ai. 

Next, construct a complete binary tree as follows, 

1. Prepare a complete binary tree with enough height 
h, this height h should satisfy 2"^-^ <N < 2''"^ 

2. Label each node with T^, which denotes the n*^ 
value at level k. The root node is labeled by T/*. 
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has the sum of all A , that is, 



FIG. 3: (Color online) Computation of the value of Ai. (a) 
Filled gray circles represent neighboring particles with radius 
a and large circles are shadows of them (we call them 'shadow 
circles') with radius 2cr. (b) An enlargement. The shaded area 
is the area into which the particle i is allowed to move. This 
figure is made by removing overlaps of shadow circles from 
the trial circle of radius as centered on the chosen particle, 
(c) The survived vertex points. The center of the trial circle 
is denoted by O. The solid circles represent survived vertex 
points, which form the area Ai. With them, we can calculate 
the value of Ai. The rectangle denotes a bounding rectangle. 
Each two adjacent survived vertex points and the center point 
O form a triangle. In this example there are five survived ver- 
tex points, and consequently five triangles to consider, (d) To 
calculate a portion of Ai, the area within each triangle formed 
by survived vertex points and O is calculated. The survived 
vertex points are the intersection points of the shadow circles 
or the trial circles, and here are denoted by Pi and P2. The 
center of the shadow particle is S. To find the shaded area 
a Monte Carlo procedure is performed in the shaded area of 
either (c) or (d), and only survived points generated in the 
shaded area are used as the new location for the new point 
O. 



N 



(7) 



Using this tree, we can choose a particle with the proba- 
bility proportional to Ai in the following way. 

1. fc ^ 1, i ^ 1. 

2. Prepare a random number r uniform on (0, T*). 
1 



3. 



2i 
2i 



otherwise 



A. k^i-1. 



5. if fc > 1 then go to|21 

Consequently, choosing the bottom node requires h — 1 
random numbers. 

After the above processes, the final value of i indicates 
the index of the particle to move. The order of this search 
algorithm is 0(log A^). When the position of particle i is 
moved, the value of Ai is also modified. We only have 
to update part of this tree for the chosen particle and 
its neighbors. The order of this update is also 0(log A^), 
which is much faster than 0{N) of the direct implemen- 
tation. Details to implement the complete binary tree 
search method are described in the appendix. 
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A node labeled T,j+^ has branches leading to two 
nodes T^^_^ and r|„. 

3. Associate every bottom node Tl with the value of 
area Ai. If the number of bottom nodes 2'*"^ is 
larger than A^, the rest of the nodes are associated 
with zero, namely, T/ = (i > A^). 

4. Associate nodes at higher levels [k > 1) recursively 
with the sum of the values associated with its two 
children, namely, T„^+i = T|„_i + T^^. 

A sample of a compete binary tree is shown in Fig. 0] 
Each node has the value and the value of each node 
at level fc -|- 1 is the sum of the values of its two children 
nodes at level k. The root node, which is Tj* in Fig. 21 



FIG. 4: A complete binary tree search. An example of A = 8 
{h = 4) is shown. The value of a node T/" is the sum of the 
values of its two child nodes, namely, Ti = T^~^ 4- T^^lX- 
After construction of the tree, we use h — 1 random numbers 
to choose a bottom node. This bottom node is associated with 
the index of the particle that will be moved in this algorithmic 
step of the RFMC method. 



E. Find a new position of the particle 

After choosing a particle, we have to choose the new 
position for it. It is difficult to choose a point uniformly 
from the points allowed to move into, since its shape is 
generally very complicated (see Fig. Therefore, we 
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have chosen to choose the new position using a Monte 
Carlo rejection method. Namely we generate a random 
position uniformly over some bounding area that includes 
all of the area Ai. [Such as the dashed circle of radius Cs 
in Figl^Ja) or the rectangle in Fig|2Ic).] If this point is 
not in Ai it is rejected and another uniformly distributed 
random point over the bounding area is generated. The 
first point not be rejected is the new position of the par- 
ticle, since it is in the area Ai, and this point is used 
as the new center for the particle. This completes one 
algorithmic step of the RFMC method. 

A value of area Ai is very small compared to the trial 
circle at high density, and hence the Monte Carlo trial to 
find the new position of the particle to be moved became 
very inefficient. To improve this, it is effective to limit the 
trial area for the Monte Carlo by making the bounding 
area very close to the area Ai. We outline two different 
survived point methods, but have only implemented the 
first. 

For the first method, the one actually implemented 
in this paper see Fig. O (c). The solid rectangle is a 
bounding rectangle which includes the area Ai. It is 
easy to obtain the bounding rectangle with the survived 
vertex points. With the set of survived vertex points 
{{xi^Ui)}, a diagonal line of the bounding rectangle is 
from (min {x;}, min {j/i}) to (max {xi}, max {i/i}). Then 
we can perform Monte Carlo trials for a new position 
within only this rectangle. The area of the rectangle is 
on the same order of Ai, so the probability of success to 
obtain the new position is drastically improved compared 
with the direct search over the trial circle. 

An alternative method is first use a random number 
to decide which of the triangles formed with point O and 
two adjacent survived vertex points the survived point 
will fall into. This is done analytically since the area of 
each of these triangles with removed shadow circle chords 
have been already calculated. Then the shortest side 
formed with point O and the two survived vertex points 
(say SP2 in Fig.|21(d)) is lengthened to be equal to the 
longest side {SPi in Fig. 121(d)). The random trial point 
is then generated within the section of the circle with 
a radius equal to the longest side {SPi in Fig. 01 (d)). 
Then the point becomes the survived point used for the 
new location of point O if the trial point is within the 
shaded area. Otherwise, this procedure repeats in the 
same extended circular section until a survived point is 
found. 



III. RESULTS 

A. Calculation of Ai 

In order to test our method to calculate Ai described in 
Sec. Ill Cl the values of Ai were also evaluated by a Monte 
Carlo sampling (Amc) with trial points uniformly drawn 
over the trial circle. The density of the system p is de- 
fined to be p = ANa'^ /L'^ with the number of particles iV, 



the radius of the particles a and the linear system size 
L, respectively. Throughout this study, the number of 
particles N is set to be 23288 and the periodic boundary 
condition is taken for both axes. The number of the gen- 
erated configurations were 3000, and 10^ MC trial points 
are taken for each of the configurations to evaluate its 
area. The density of the system is fixed at p = 0.9. The 
result is shown in Fig. |S1 The area Ai is normalized by 
the area of the trial circles (see Fig. jJJ. The difference 
between the MC and our survived point method is less 
than 0.01% for all areas, which is the same order as the 
statistical error of our MC method. The standard devia- 
tion of the MC area calculation is determined by dividing 
the data into 10 groups, each including 10^ samples. This 
result shows that the value of Ai is properly calculated 
by our method. 
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FIG. 5: Comparison of calculated Ai between our survived 
points method and that calculated by a more straightforward 
MC method. Units on both axes are dimensionless. The cir- 
cles denote the ratio of j4]^p]y[Q to the Ay^Q. The number 
of configurations is 3000 at p = 0.9 and 10® independent sam- 
ples are averaged for each configuration. The solid lines are 
the standard deviation obtained from the jack-knife procedure 
described in the text. 



B. Time evolution 

To compare the dynamics between the standard MC 
and the RFMC, we observed the time evolution of the 
six- fold bond-orientational order parameter ,2Qj . The 
parameter is defined to be 

<^6 = (exp(6i0)) , (8) 

with the bond angle 9 which has a definition described 
in Fig. (HI The average is taken for all pairs of neigh- 
boring particles. The parameter becomes 1 when all 
particles are located on the points of a hexagonal grid, 
and it becomes when the particle location is completely 
disordered. Therefore (f)e describes how close the system 
is to the perfect hexagonal packing. The neighbors in 
an off-lattice model are strictly defined with the Voronoi 
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construction [Sj , which is a very time-consuming method. 
In this paper, two particles separated by a distance less 
than 2.6(7 are defined as neighbors. We confirmed that 
the value of (/)g is approximately the same value as the 
value obtained with the Voronoi construction. At the 




1 

0.9 
0.8 
0.7 

s 

0.5 
0.4 
0.3 
0.2 
0.1 

















p=0.60 




0.70 




0.80 




0.90 




1.00 





20 



40 



60 



80 



100 



Area to search 



Definition of the angle 



FIG. 6: The definition of the neighboring particles and bond 
angle 6. Two particles separated by a distance less than R 
are defined as neighbors. Here, R is set to be 2.6a with the 
radius a of particles. The bond angle 9 is defined to be an 
angle between the bond connecting neighboring particles with 
respect to an arbitrary, but fixed global axis. 



FIG. 7: Relaxation behaviors of the bond-orientational or- 
der. Solid lines are the results of the standard MC and 
symbols (circles, triangles and squares) are the results of the 
RFMC. The data intervals between accepted updates for the 
RFMC algorithm becomes longer at high density, while the 
data keeps the behavior of the standard MC algorithm. 



beginning of the simulation, the particles are set up in a 
perfect hexagonal order, namely, (jjeit = 0,p) = 1. The 
order parameter (j)Q starts to relax to the value of the 
equilibrium state. With this nonequilibrium relaxation 
(NER) behavior of order parameters, critical points and 
critical exponents of various phase transitions can be de- 
termined accurately |2lj|, |^ H^]. This method is 
called a NER method. Watanabe et al. 0, ^ stud- 
ied two-dimensional melting based on the NER method 
for the Kostcrlitz-Thouless transition [2^ by observing 
the relaxation behavior of <j)Q. Therefore, the following 
time evolutions of 06 contains information about the two- 
dimensional melting transition. 

Time evolutions of (pe are shown in Fig. |7| Solid 
lines are results of the standard MC simulation and sym- 
bols (circles, triangles and squares) are the results of the 
RFMC. Fig. [3 shows that both behaviors are equivalent 
for the two methods. This is essentially a check of the 
program implementation, since the physical dynamic is 
the same for both the MC and the RFMC methods. 



C. Efficiency 

The computational times required to achieve 1000 ac- 
cepted MC steps are shown in Fig.|SIa). Configurations 
are started from the perfect hexagonal configuration and 
both measurements are started after 10*^ MC steps. All 
simulations are performed on an Intel Xeon 2.4 GHz com- 
puter. While the computational time of the RFMC (open 
circles) is almost constant, a longer computational time 
is required for the standard Monte Carlo (solid circles) at 
higher density. It shows that the RFMC is more efficient 
at high densities, in spite of the additional bookkeeping 
involved in the RFMC method (so one algorithmic step 
takes much longer than one standard MC step). 



The CPU-time ratio of the standard MC to the RFMC 
is shown in Fig. IHl^b) . The data are shown as a function 
of 1/e, where e = (pcp — p)/ Pc the closest density is pcp 
and the density of the system is p. The CPU-time ratio, 
which is the eflaciency of the RFMC compared to that of 
the standard MC, diverges as e~^. 



IV. SUMMARY AND DISCUSSION 

We constructed a rejection-free Monte Carlo algorithm 
for the hard-disk system. This method conserves the 
property of the dynamic behavior of the original Monte 
Carlo method. In other words, the time scales will 
not depend on the density, but are rather set by some 
Brownian-motion type of dynamic for all densities. An 
estimate of the time scales between the MC and physical 
time can thus be obtained by setting the mean-free path 
of an isolated particle to be proportional to the value as ■ 
Note that strictly this is only true in the limit CTs — > 0, 
but it should be a reasonable approximation for a small 
finite Us- 

We also find that for a fixed value of CTs, the RFMC 
method is more efficient at high density. Therefore, 
the RFMC method should be useful for studies of two- 
dimensional solids or studies of high-density glass mate- 
rials. It may also be possible to make the algorithm even 
more efficient by utilizing the ideas of absorbing Markov 
chains (for the MCAMC method for discrete state spaces 
see [3 and references therein). Increased algorithmic 
efficiencies for the Monte Carlo dynamics of hard disks 
could be useful to further test physical phenomena using 
hard disk systems, such as for example the relationship 
between fiuctuations and dissipation of work in a Joule 
experiment '2^1 . 
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Appendix 

The complete binary tree search can be implemented 
with an one- dimensional array. To make it simple, let 
the number of particles N be 2^~^. The tree with height 
h requires an array a{i) with size 2A'' — 1. First, associate 



a(l) a(2) o(3) o(4) a(5) q(6) a(7) 



a.(N) 



a(2N - 1) 



each bottom node with a corresponding value as 

a{N + i-l)^Ai (i = 1,2, •••,7V), (9) 

which corresponds to T/ «— Aj. Next, associate parent 
nodes recursively as 

i ^ N 
While i ^ 

a{i) ^ a{2i) + a{2i+ 1) 
i <— i — 1 
next i 

Using this array, we can pick up particle i with the prob- 
ability proportional to A, as, 

i ^ 1 

While i < N 

Prepare a uniform random number r on (0, a{i)) 



i if r < a{2i) 

« <— 2i + 1 otherwise 



next i 

i ^ i - N + 1. 

After the above procedure, we obtain the index i of the 
chosen particle. When the value of Ai is changed, the tree 
should be updated. The update process is as follows, 

a{N + i-l)^Ai 

i + N 
2 

While i 7^ 1 

ail) ^ a{2i) + a{2i+ 1) 
i ^ [i/2\ 
next i. 

Note that, when the chosen particle is moved, the ac- 
ceptance probabilities of the neighboring particles of the 
moved particle are also changed. Therefore, we have to 
perform the above process for all neighboring particles. 



rJi-i 

'2 



rpll-2 rpll-2 



rh-2 



rph- 
J-4 



Tl 



J.1 



FIG. 9: Implementation of the complete binary tree search 
with an array. The required size of the array to implement 
the tree with height h is 2^ — 1. The height h should satisfy 
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